ending = '_US';
addpath('.\functionsRA')
load('.\output\results_est_1_US.mat')
load('.\input\dists_US.mat','ny','years')

% graph with comparison of the models
models = 1;                                                                 % here is a vector with the models that will be compared in the graph
qs = zeros(12*ny,size(models,1));
for i = 1:size(models,1)
    m = models(i);
    qs(:,i) = reshape(Qs(:,:,m),ny*12,1);
end
qs = qs(qs>0); lth = size(qs,1);                                            % we stop in Mar 2021
p = plot((1:lth),100*qs(:,1),(1:lth),100*mean(qs).*ones(lth,1),'--','LineWidth',2);
p(1).Color = '#A2142F'; %p(2).Color = '#0072BD'; p(3).Color = '#77AC30'; p(4).Color = '#7E2F8E';
p(2).Color = '#A2142F'; %p(6).Color = '#0072BD'; p(7).Color = '#77AC30'; p(8).Color = '#7E2F8E';
title('Result of Minimisations for US (# parameters)')
legend('Model 1 (k=6/month)','Location','northwest')
axis([1 lth 0 10])
set(gca,'Xtick',(1:12:ny*12),'Ytick',(0:1:10),'XtickLabel',years,'XTickLabelRotation',45);
ax = gca; ax.FontSize = 16; grid on
ax.YAxis.TickLabelFormat = '%g%%';
set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
print('-bestfit','.\figures matlab\comparison_US.pdf','-dpdf')

mdl = 1;
pis_graph = reshape(pis{mdl,1},6,12*ny)'; pis_graph(lth+1:end,:) = [];
name = ['model_' num2str(mdl) '' ending '\model_' num2str(mdl)];
% graph with probabilities for the A matrix
p2 = plot((1:lth),100*pis_graph,'LineWidth',2);
title(['Probabilities for US, Model (' num2str(mdl) ')'])
if mdl == 1
    legend('\pi_{L}','\pi_{H}','\pi_{nL}','\pi_{nH}','\pi_{nn}','\pi_{mr}','Location','northwest')    
end
axis([1 lth 0 85])
set(gca,'Xtick',(1:12:ny*12),'Ytick',(0:10:80),'XtickLabel',years,'XTickLabelRotation',45);
ax = gca; ax.FontSize = 16; grid on
ax.YAxis.TickLabelFormat = '%g%%';
set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
print('-bestfit',['.\figures matlab\' name '' ending '_probs.pdf'],'-dpdf')

% graphs with the marginal densities in the model and in the data
xgrid = (1:8)';
xlabels = {'<-1%','-1-0%','0-1%','1-2%','2-3%','3-4%','4-5%','>5%'};
fs_d_graph = reshape(fs_data,8,12*ny)'; fs_d_graph(lth+1:end,:) = []; 
g5s_d_graph = reshape(squeeze(gs_data(:,1,:,:)),8,12*ny)'; g5s_d_graph(lth+1:end,:) = [];
g10s_d_graph = reshape(squeeze(gs_data(:,2,:,:)),8,12*ny)'; g10s_d_graph(lth+1:end,:) = [];
fs_m_graph = reshape(fs_model{mdl,1},8,12*ny)'; fs_m_graph(lth+1:end,:) = []; 
g5s_m_graph = reshape(squeeze(gs_model{mdl,1}(:,1,:,:)),8,12*ny)'; g5s_m_graph(lth+1:end,:) = [];
g10s_m_graph = reshape(squeeze(gs_model{mdl,1}(:,2,:,:)),8,12*ny)'; g10s_m_graph(lth+1:end,:) = [];

name = ['model_' num2str(mdl) '' ending '\marginals\model_' num2str(mdl)];
months = {'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'}';
ms = repmat((1:12),1,ny)'; ys = repmat((1:ny),12,1); ys = ys(:);
for i = 1:lth
plot(xgrid,100*[fs_d_graph(i,:); fs_m_graph(i,:)],'LineWidth',2)
title(['Avg of Margs t+6 to t+10 in ' months{ms(i)} '/' years{ys(i)} ', US, Model (' num2str(mdl) ')'])
legend('Data','Model','Location','northwest')
axis([1 xgrid(end) 0 50])
set(gca,'Xtick',xgrid,'Ytick',(0:5:50),'XtickLabel',xlabels);
ax = gca; ax.FontSize = 16; grid on
ax.YAxis.TickLabelFormat = '%g%%';
set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
print('-bestfit',['.\figures matlab\' name '' ending '_marginals_' years{ys(i)} months{ms(i)} '.pdf'],'-dpdf')
end

% graphs with the ZC t+5 probabilities in the model and in the data
name = ['model_' num2str(mdl) '' ending '\avg5\model_' num2str(mdl)];
for i = 1:lth
plot(xgrid,100*[g5s_d_graph(i,:); g5s_m_graph(i,:)],'LineWidth',2)
title(['Dist of Avg Inf t+5 in ' months{ms(i)} '/' years{ys(i)} ', US, Model (' num2str(mdl) ')'])
legend('Data','Model','Location','northwest')
axis([1 xgrid(end) 0 50])
set(gca,'Xtick',xgrid,'Ytick',(0:5:50),'XtickLabel',xlabels);
ax = gca; ax.FontSize = 16; grid on
ax.YAxis.TickLabelFormat = '%g%%';
set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
print('-bestfit',['.\figures matlab\' name '' ending '_avg5_' years{ys(i)} months{ms(i)} '.pdf'],'-dpdf')
end

% graphs with the ZC t+10 probabilities in the model and in the data
name = ['model_' num2str(mdl) '' ending '\avg10\model_' num2str(mdl)];
for i = 1:lth
plot(xgrid,100*[g10s_d_graph(i,:); g10s_m_graph(i,:)],'LineWidth',2)
title(['Dist of Avg Inf t+10 in ' months{ms(i)} '/' years{ys(i)} ', US, Model (' num2str(mdl) ')'])
legend('Data','Model','Location','northwest')
axis([1 xgrid(end) 0 50])
set(gca,'Xtick',xgrid,'Ytick',(0:5:50),'XtickLabel',xlabels);
ax = gca; ax.FontSize = 16; grid on
ax.YAxis.TickLabelFormat = '%g%%';
set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
print('-bestfit',['.\figures matlab\' name '' ending '_avg10_' years{ys(i)} months{ms(i)} '.pdf'],'-dpdf')
end

% graph with the prob of average inflation from t+6~t+10 being a disaster
name = ['model_' num2str(mdl) '' ending '\model_' num2str(mdl)];
pi_lim = [(-1:5)'; Inf];
pi_bar = [-2; (-0.5:4.5)'; 6];
prob_avg_dis = zeros(lth,2); prob_avg_dis04 = prob_avg_dis;
load('.\input\dists_US.mat','infl')
for i = 1:lth
    bin = (pi_lim>=infl(i)); bt = find(bin,1);
    A = Amatrix(pis_graph(i,:)',mdl);
    g = avg_inf_6to10(A,bt,pi_bar,pi_lim);
    prob_avg_dis(i,:) = [g(1) g(end)];

    prob_avg_dis04(i,:) = [g(1)+g(2) g(end-1) + g(end)];
end
plot((1:lth),100*prob_avg_dis,'LineWidth',2)
title(['Prob of Avg Infl t+6~10 is a disaster, US, Model(' num2str(mdl) ')'])
legend('Deflation','Hyperinflation','Location','northeast')
axis([1 lth 0 20])
set(gca,'Xtick',(1:12:ny*12),'Ytick',(0:2:20),'XtickLabel',years,'XTickLabelRotation',45);
ax = gca; ax.FontSize = 16; grid on
ax.YAxis.TickLabelFormat = '%g%%';
set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
print('-bestfit',['.\figures matlab\' name '' ending '_prob_avgdis.pdf'],'-dpdf')

% graph with the prob of average inflation from t+6~t+10 being a disaster with thresholds of 0 and 4 instead of -1 and 5
plot((1:lth),100*prob_avg_dis04,'LineWidth',2)
title(['Prob of Avg Infl t+6~10 is a disaster, US, Model(' num2str(mdl) ')'])
legend('Deflation','Hyperinflation','Location','northeast')
axis([1 lth 0 31])
set(gca,'Xtick',(1:12:ny*12),'Ytick',(0:2:31),'XtickLabel',years,'XTickLabelRotation',45);
ax = gca; ax.FontSize = 16; grid on
ax.YAxis.TickLabelFormat = '%g%%';
set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
print('-bestfit',['.\figures matlab\' name '' ending '_prob_avgdis04.pdf'],'-dpdf')

R2(R2(:,1)==0,:) = [];
% graph on R2 and RMSE
p = plot((1:lth),100*R2(:,1),'-',(1:lth),R2(:,1)*0+100*mean(R2(:,1)),'r--','LineWidth',2);
p(1).Color = '#A2142F';
p(2).Color = '#A2142F';
title(['Model Fit, US, Model(' num2str(mdl) ')'])
axis([1 lth 0 100])
set(gca,'Xtick',(1:12:ny*12),'Ytick',(0:10:100),'XtickLabel',years,'XTickLabelRotation',45);
ax = gca; ax.FontSize = 16; grid on
ax.YAxis.TickLabelFormat = '%g%%';

yyaxis right
p = plot((1:lth),100*qs,'-',(1:lth),qs*0+100*mean(qs),'--','LineWidth',2);
p(1).Color = '#0072BD';
p(2).Color = '#0072BD';
axis([1 lth 0 20]); ax = gca; ax.FontSize = 16; box on; grid on
ax.YAxis(2).Color = [0 0 0]; ax.YAxis(2).TickLabelFormat = '%g%%';
set(gcf,'Units','centimeters','PaperSize',[17 12])
set(gca,'Ytick',(0:2:20),'Units','centimeters','position',[1.8 3 11.3 7.3]); % [left bottom width height])
legend('R2, lhs','mean R2','RMSE, rhs','mean RMSE','Location','southoutside','NumColumns',4,'Units','centimeters','position',[1.45 0.3 12 0.7])
print('-bestfit',['.\figures matlab\' name '' ending '_fit.pdf'],'-dpdf')
